%clear all; clc;close all;
%Updated: 11/30/2021

%plot policy functions and value functions
a   =load('agrid2.txt');
b   =load('bgrid2.txt');

vcon0=load('gvcon.txt');
vaut0=load('gvaut.txt');

p0  =load('gp.txt');
h0  =load('gh.txt');
ibp0=load('gibp.txt');
ct0 =load('gct.txt');
ca0 =load('gca.txt');
gn0 =load('ggn.txt');
tau0=load('gtau.txt');
w0  =load('gw.txt');
def0=load('gdef1.txt');
q0  =load('gq.txt');

paut0  =load('gpaut.txt');
haut0  =load('ghaut.txt');
ctaut0 =load('gctaut.txt');
gnaut0 =load('ggnaut.txt');
tauaut0=load('gtauaut.txt');
waut0  =load('gwaut.txt');

fm0     =load('fm_ss.txt');
rchange0=load('rchange_ss.txt');
derivb0=load('derivativeb_ss.txt');
bpvalue30=load('bpoptimal_ss.txt');
bpvalue20=load('bpgn2_ss.txt');
qvalue30=load('qoptimal_ss.txt');
qvalue20=load('qgn2_ss.txt');

na=rows(a);
nk=1;
ns=na/nk;
nb=rows(b);

alpha  = 0.75;
pm     = 1;
gamma  = 0.43;
omegac = 0.3;
mu     = 1;
lambda = 0.184;
r      = 0.02;

atrans = ones(nb,1)*a';

for i=1:na;
    for j=1:nb;
        vcon(i,j)  =vcon0((i-1)*nb + j);
        
        p(i,j)  =p0((i-1)*nb + j);
        h(i,j)  =h0((i-1)*nb + j);
        ibp(i,j)=ibp0((i-1)*nb + j);
        ct(i,j) =ct0((i-1)*nb + j);
        ca(i,j) =ca0((i-1)*nb + j);
        gn(i,j) =gn0((i-1)*nb + j);       
        tau(i,j)=tau0((i-1)*nb + j);       
        w(i,j)  =w0((i-1)*nb + j);   
        def(i,j)=def0((i-1)*nb + j);   
        q(i,j)  =q0((i-1)*nb + j);   
        rec(i,j)=w(i,j)*h(i,j)*tau(i,j);
    end;
end;

bp = b(ibp);

for i=1:na;
    vaut(i)  =vaut0(i);
    
    paut(i)  =paut0(i);
    haut(i)  =haut0(i);
    ctaut(i) =ctaut0(i);
    gnaut(i) =gnaut0(i);       
    tauaut(i)=tauaut0(i);       
    waut(i)  =waut0(i);   
    recaut(i)=waut(i)*haut(i)*tauaut(i);
	wwaut(i) = alpha*paut(i).*haut(i).^(alpha-1); 
end;

cn       = h.^alpha-gn;
cnaut    = haut.^alpha-gnaut;
const    = (omegac.*ct.^(-mu)+(1-omegac).*cn.^(-mu)).^(-1/mu);
constaut = (omegac.*ctaut.^(-mu)+(1-omegac).*cnaut.^(-mu)).^(-1/mu);

for i=1:na;
    temp1 = find(ct(i,:) > 0);
    temp2 = find(cn(i,:) > 0);
    temp3 = find(p(i,:) > 0);

    iminct(i) = temp1(1);
    imincn(i) = temp2(1);
    iminp(i) = temp3(1);
    
    imin(i)  = max(iminct(i),max(imincn(i),iminp(i)));
    
    v(i,1:imin(i)-1)  = NaN; 
    p(i,1:imin(i)-1)  = NaN; 
    bp(i,1:imin(i)-1) = NaN; 
    ibp(i,1:imin(i)-1)= NaN; 
    ctE(i,1:imin(i)-1) = NaN; 
    ctU(i,1:imin(i)-1) = NaN; 
    cnE(i,1:imin(i)-1) = NaN; 
    cnU(i,1:imin(i)-1) = NaN; 
    constE(i,1:imin(i)-1) = NaN;     
    constU(i,1:imin(i)-1) = NaN;     
    gn(i,1:imin(i)-1) = NaN; 
    h(i,1:imin(i)-1)  = NaN; 
    ca(i,1:imin(i)-1)  = NaN; 
    tau(i,1:imin(i)-1)= NaN; 
	rec(i,1:imin(i)-1)= NaN;
	ww(i,1:imin(i)-1)= NaN;
	q(i,1:imin(i)-1)  = NaN;
  
end;

for i=1:ns;
for j=1:nk;        
    vcon2(i,j,:)  = vcon((i-1)*nk+j,:); 
    p2(i,j,:)  = p((i-1)*nk+j,:);  
    bp2(i,j,:) = bp((i-1)*nk+j,:);  
    ibp2(i,j,:)= ibp((i-1)*nk+j,:);  
    ct2(i,j,:) = ct((i-1)*nk+j,:);  
    cn2(i,j,:) = cn((i-1)*nk+j,:);  
    const2(i,j,:) = const((i-1)*nk+j,:);         
    gn2(i,j,:) = gn((i-1)*nk+j,:);  
	ca2(i,j,:) = ca((i-1)*nk+j,:);  
    h2(i,j,:)  = h((i-1)*nk+j,:);   
    tau2(i,j,:)= tau((i-1)*nk+j,:);  
	rec2(i,j,:)= rec((i-1)*nk+j,:); 
    def2(i,j,:)= def((i-1)*nk+j,:);  
	q2(i,j,:)  = q((i-1)*nk+j,:); 
	ww2(i,j,:) = alpha*p2(i,j,:).*h2(i,j,:).^(alpha-1); 
    

    vaut2(i,j)  = vaut((i-1)*nk+j); 
    paut2(i,j)  = paut((i-1)*nk+j);  
    ctaut2(i,j) = ctaut((i-1)*nk+j);  
    cnaut2(i,j) = cnaut((i-1)*nk+j);  
    constaut2(i,j) = constaut((i-1)*nk+j);          
    gnaut2(i,j) = gnaut((i-1)*nk+j);  
    haut2(i,j)  = haut((i-1)*nk+j);  
    tauaut2(i,j)= tauaut((i-1)*nk+j);  
	wwaut2(i,j) = alpha*paut2(i,j).*haut2(i,j).^(alpha-1); 
end;
end;

nettax2    = p2.*gn2-ca2;           %net taxes
nettaxaut2 = paut2.*gnaut2;         %net taxes

fiscal2   = tau2 - p2.*gn2;           %fiscal surplus
fiscalaut2=tauaut2-paut2.*gnaut2;     %fiscal surplus

nfiscal2   = nettax2 - p2.*gn2;           %fiscal surplus net of transfers
nfiscalaut2=nettaxaut2-paut2.*gnaut2;     %fiscal surplus net of transfers

lowa = round(2*ns/6);
meda = round(ns/2);
hgha = round(4*ns/6);

b =-b/(lambda+r);  %change of sign
a = exp(a);

%Compute fiscal multipliers an borrowing costs
for i=1:nb;
for j=1:na;
	fm(j,i)      =fm0((i-1)*na + j);     
  	rchange(j,i) =rchange0((i-1)*na + j);     
  	derivb(j,i)  =-derivb0((i-1)*na + j);     
  	bpvalue2(j,i)=bpvalue20((i-1)*na + j);     
  	bpvalue3(j,i)=bpvalue30((i-1)*na + j);     
  	qvalue2(j,i)=qvalue20((i-1)*na + j); 
    qvalue2(j,i)=100*((1+qvalue2(j,i)^(-1)-lambda)/(1+r)-1);
  	qvalue3(j,i)=qvalue30((i-1)*na + j);   
    qvalue3(j,i)=100*((1+qvalue3(j,i)^(-1)-lambda)/(1+r)-1);
	difr(j,i)= -((1+qvalue3(j,i)^(-1)-lambda)/(1+r)-1) ;
    difr(j,i)=bpvalue2(j,i)-bpvalue3(j,i);
    if (def(j,i)==0);fm(j,i)=NaN;rchange(j,i)=NaN;derivb(j,i)=NaN;bpvalue2(j,i)=NaN;...
            bpvalue3(j,i)=NaN;qvalue2(j,i)=NaN;qvalue3(j,i)=NaN;difr(j,i)=NaN;end;     
end;
end;

if strcmp(label_graph,'fiscal_multipliers');

%income index =11 (mean income)

%debt indices
%default IR: 
%=122 mean debt,        value grid = -0.2261
%=196 0.5*mean debt,    value grid = -0.1136
%=159 0.75*mean debt,   value grid = -0.1698

%extra info: bss = 108, value grid = -0.2474

%risk-free debt IR:
%%=1888 mean debt,      value grid = -0.226
%=1944 0.5*mean debt,   value grid = -0.114
%=1916 0.75*mean debt,  value grid = -0.17
%=1774 2*meandebt       value grid = -0.454
%=1877 1.1*mean debt    value grid = -0.248
%=1973 0.25*mean debt,  value grid = -0.056
%=2001 zero,            value grid = 0


%debt indices for fiscal mutipliers:
% 1.1*mean debt  == 107       value grid = -0.2489
% mean debt      == 122       value grid = -0.2261
% 0.75*mean debt == 159       value grid = -0.1698 

hghb = 107;     
xmin = 0.95;
xdif = 0.05;
xmax = 1.15;
xf   = 16; %fontsize

[f1,f2]=find(fm(:,hghb)>1);
    fm1=zeros(size(fm(:,hghb)));
    fm1(:,:)= NaN;
    fm1(f1)=fm(f1,hghb);
    fm2= fm(:,hghb);
    fm2(f1)=NaN;
    
NameArray1 = {'LineStyle','Color','LineWidth'};
ValueArray1 = {'-';[0.078  0.1686  0.549];2}';

nameFig = char([char(label_model) '_' char(label_graph) '_yt_paper2']);
figure('name',nameFig);
P=plot(a,fm1,'-b',a,fm2,'-b','LineWidth',2);
set(P,NameArray1,ValueArray1);grid on;
AX1 = gca;box on;
set(AX1,'YTick',-1:1:3,'FontSize',xf);
axis(AX1, [xmin xmax -1.1 3]);
%title('fiscal mutipliers');
set(AX1,'XTick',xmin:xdif:xmax,'FontSize',xf,'XLim',[xmin xmax]);
xlabel('y^T','FontSize',xf);
set(gcf, 'Units', 'Inches', 'Position', [2, 2, 7, 5]); %[initial x position, initial y position, width, height]
saveas(gcf, fullfile(plot_path, nameFig),format_chart);

nameFig = char([char(label_model) '_' char(label_graph) '_bc_paper2']);
figure('name',nameFig);
P=plot(a,derivb(:,hghb),'-b','LineWidth',2);
set(P,NameArray1,ValueArray1);grid on;
AX1 = gca;box on;
set(AX1,'YTick',0:0.5:1,'FontSize',xf);
set(AX1,'XTick',xmin:xdif:xmax,'FontSize',xf,'XLim',[xmin xmax]);
axis(AX1, [xmin xmax -0.05 1]);
xlabel('y^T','FontSize',xf);
set(gcf, 'Units', 'Inches', 'Position', [2, 2, 7, 5]); %[initial x position, initial y position, width, height]
saveas(gcf, fullfile(plot_path, nameFig),format_chart);

return
end;

%put together policy fcns of repayment and autarky
mat_ones = ones(size(def));
risk_free_debt=0;

if (risk_free_debt==1);
    p3  =p;
    h3  =h;
    ct3=ct;
    cn3=cn;
    ca3 =ca;
    gn3 =gn;       
    tau3=tau;       
    w3  =w;   
    rec3=rec;
else;
    p3  =def.*p+(mat_ones-def).*(paut'*ones(1,nb));
    h3  =def.*h+(mat_ones-def).*(haut'*ones(1,nb));
    ct3=def.*ct+(mat_ones-def).*(ctaut'*ones(1,nb));
    cn3=def.*cn+(mat_ones-def).*(cnaut'*ones(1,nb));
    ca3 =def.*ca+(mat_ones-def).*0;
    gn3 =def.*gn+(mat_ones-def).*(gnaut'*ones(1,nb));       
    tau3=def.*tau+(mat_ones-def).*(tauaut'*ones(1,nb));       
    w3  =def.*w+(mat_ones-def).*(waut'*ones(1,nb));   
    rec3=def.*rec+(mat_ones-def).*(recaut'*ones(1,nb));
end;

[i1,i2] = find(def>0);
q3  =NaN; bp3 =NaN; sp3 =NaN; gn3 = NaN; h3=NaN; p3=NaN;

for i=1:na;
for j=1:nb;
    if (def(i,j)>0); q3(i,j)=q(i,j);bp3(i,j)=bp(i,j);sp3(i,j)=100*((1+q3(i,j).^(-1)-lambda)./(1+r)-1);p3=p;h3=h;gn3=gn;tau3=tau;w3=w;rec3=rec;
    else
        q3(i,j)=NaN;bp3(i,j)=NaN;sp3(i,j)=NaN;p3=NaN;h3=NaN;gn3=NaN;tau3=NaN;w3=NaN;rec3=NaN;
    end;
end;
end;

nettax3    = p3.*gn3-ca3;           %net taxes
fiscal3    = tau3 - p3.*gn3;        %fiscal surplus
nfiscal3   = nettax3 - p3.*gn3;     %fiscal surplus net of transfers

%value functions: repayment vs autarky
vcon21(:,:) = vcon2(:,1,:);
vcon22(:,:) = vcon2(:,1,:);
vaut21(:,:) = vaut2(:,1);
vaut22(:,:) = vaut2(:,1);

figure;
plot(b(1:nb),vcon21(lowa,:),'-r',b(1:nb),vcon21(meda,:),':r',b(1:nb),...
    vaut21(lowa,:)*ones(size(b)),'-b',b(1:nb),vaut21(meda)*ones(size(b)),':b','LineWidth',1.5);
xlabel('b');title('value functions');
ylim([-16 -13]);legend('rep low y^T','rep high y^T','aut low y^T','aut high y^T')

%% Plot variables as function of yT for two debt levels
%optimal policies and prices: repayment vs autarky

if strcmp(label_graph,'policy_fcns');

NameArray = {'LineStyle','Color','LineWidth'};
NameArray0 = {'LineStyle','Color','LineWidth'};
ValueArray0 = {'--','-';[1  0.4  0.4],[0.078  0.1686  0.549];1.5,1.5}';
ValueArray = {'--','-';[1  0.4  0.4],[0.078  0.1686  0.549];1.5,1.5}';
ValueArray2 = {'--','-','--','-';[1  0.4  0.4],[0.078  0.1686  0.549],[1  0.4  0.4],[0.078  0.1686  0.549];1.5,1.5,1.5,1.5}';
ValueArray4 = {'--','--','-','-';'red','blue','red','blue';1.5,1.5,1.5,1.5}';
ValueArray5 = {':',':';[1  0.4  0.4],[0.078  0.1686  0.549];1.5,1.5}';

endb = 2;
lowb = 159;
hghb = 122;

nameFig = char([char(label_model) '_' char(label_graph) '_yt']);
figure('name',nameFig);
subplot(5,2,1);P = plot(a,100*squeeze(1-h3(:,lowb)),a,100*squeeze(1-h3(:,hghb)));
title('unemployment','Fontweight','normal');
AX1 = gca;box on;ylabel('%');
set(AX1,'YTick',0:20:40);
axis(AX1, [a(1) a(end) -2 40]);
set(P,NameArray,ValueArray);
subplot(5,2,2);P = plot(a,squeeze(p3(:,lowb)),a,squeeze(p3(:,hghb)));
title('p^{N}','Fontweight','normal'); 
%lgd=legend('low y^T','high y^T');lgd.Location='SouthEast';
legend('low debt','high debt','Location','SouthEast','Orientation','vertical');
AX1 = gca;box on;
set(AX1,'YTick',3.5:0.5:4.5);
axis(AX1, [a(1) a(end) 3.5 4.5]);
set(P,NameArray,ValueArray);
subplot(5,2,3);P = plot(a,squeeze(gn3(:,lowb)),a,squeeze(gn3(:,hghb)));
title('g^{N}','Fontweight','normal'); 
AX1 = gca;box on;
set(AX1,'YTick',0.1:0.1:0.3);
axis(AX1, [a(1) a(end) 0.1 0.3]);
set(P,NameArray,ValueArray);
subplot(5,2,4);P = plot(a,-squeeze(bp3(:,lowb))/(lambda+r),a,-squeeze(bp3(:,hghb))/(lambda+r));
title('debt','Fontweight','normal');
AX1 = gca;box on;
set(AX1,'YTick',0.8:0.2:1.2);
axis(AX1, [a(1) a(end) 0.8 1.2]);
set(P,NameArray0,ValueArray0);
subplot(5,2,5);P = plot(a,squeeze(cn3(:,lowb)),a,squeeze(cn3(:,hghb)));
title('c^{N}','Fontweight','normal');
AX1 = gca;box on;
set(AX1,'YTick',0.5:0.2:0.9);
axis(AX1, [a(1) a(end) 0.5 0.9]);
set(P,NameArray,ValueArray);
subplot(5,2,6);P = plot(a,squeeze(ct3(:,lowb)),a,squeeze(ct3(:,hghb)));
title('c^{T}','Fontweight','normal');
AX1 = gca;box on;
set(AX1,'YTick',0.7:0.2:1.1);
axis(AX1, [a(1) a(end) 0.7 1.1]);
set(P,NameArray,ValueArray);
subplot(5,2,7);P = plot(a,squeeze(w3(:,lowb)),a,squeeze(w3(:,hghb)));
title('wages','Fontweight','normal');
AX1 = gca;box on;
set(AX1,'YTick',2.5:0.5:3.5);
axis(AX1, [a(1) a(end) 2.5 3.5]);
set(P,NameArray,ValueArray);
%version 1: with bond spreads
subplot(5,2,8);P = plot(a,squeeze(sp3(:,lowb)),a,squeeze(sp3(:,hghb)));
title('spreads','Fontweight','normal');
set(P,NameArray0,ValueArray0);ylabel('%');
AX1 = gca;box on;
set(AX1,'YTick',0:4:8);
axis(AX1, [a(1) a(end) -0.5 8]);
%version 2: with current account
subplot(5,2,9);P = plot(a,-squeeze(ca3(:,lowb)),a,-squeeze(ca3(:,hghb)));
title('CA balance','Fontweight','normal');
AX1 = gca;box on;
set(AX1,'YTick',-0.1:0.1:0.2);
axis(AX1, [a(1) a(end) -0.1 0.2]);
set(P,NameArray0,ValueArray0);
%version 3: with (net) taxes
subplot(5,2,10);P = plot(a,squeeze(tau3(:,lowb)),a,squeeze(tau3(:,hghb)));
set(P,NameArray,ValueArray);
hold on;
PP = plot(a,squeeze(nettax3(:,lowb)),a,squeeze(nettax3(:,hghb)));
title('(net) taxes','Fontweight','normal');
AX1 = gca;box on;
set(AX1,'YTick',0.6:0.3:1.2);
axis(AX1, [a(1) a(end) 0.6 1.2]);
set(PP,NameArray,ValueArray5);
saveas(gcf, fullfile(plot_path, nameFig),format_chart);

%% Plot variables as function of debt for two yT levels
hr2  = h2;
cnr2 = cn2;
gnr2 = gn2;
ctr2 = ct2;
pr2  = p2;
taur2= tau2;
nettaxr2= nettax2;
fiscalr2  = fiscal2;
nfiscalr2 = nfiscal2;
bpr2 = bp2;
wwr2 = ww2;
qr2  = q2;
spr2 = 100*((1+qr2.^(-1)-lambda)./(1+r)-1);
car2 = ca2;

%optimal policies and prices: repayment vs autarky
NameArray = {'LineStyle','Color','LineWidth'};
ValueArray = {'--','--','-','-';[1  0.4  0.4],[1  0.4  0.4],[0.078  0.1686  0.549],[0.078  0.1686  0.549];1.5,1.5,1.5,1.5}';
ValueArray2 = {'--','-','--','-';[1  0.4  0.4],[0.078  0.1686  0.549],[1  0.4  0.4],[0.078  0.1686  0.549];1.5,1.5,1.5,1.5}';
ValueArray4 = {'--','--','-','-';'red','blue','red','blue';1.5,1.5,1.5,1.5}';
ValueArray5 = {':',':',':',':';[1  0.4  0.4],[1  0.4  0.4],[0.078  0.1686  0.549],[0.078  0.1686  0.549];1.5,1.5,1.5,1.5}';

NameArray0 = {'LineStyle','Color','LineWidth'};
ValueArray0 = {'--','-';[1  0.4  0.4],[0.078  0.1686  0.549];1.5,1.5}';

[ilowa_aut]=find(def2(lowa,1,:)==0);
[ilowa_rep]=find(def2(lowa,1,:)==1);
hr2(lowa,1,ilowa_aut) = NaN;
ha2(lowa,1,1:nb) = haut2(lowa,1);
ha2(lowa,1,ilowa_rep) = NaN;

[ihgha_aut]=find(def2(hgha,1,:)==0);
[ihgha_rep]=find(def2(hgha,1,:)==1);
hr2(hgha,1,ihgha_aut) = NaN;
ha2(hgha,1,1:nb) = haut2(hgha,1);
ha2(hgha,1,ihgha_rep) = NaN;

pr2(lowa,1,ilowa_aut) = NaN;
pa2(lowa,1,1:nb) = paut2(lowa,1);
pa2(lowa,1,ilowa_rep) = NaN;

pr2(hgha,1,ihgha_aut) = NaN;
pa2(hgha,1,1:nb) = paut2(hgha,1);
pa2(hgha,1,ihgha_rep) = NaN;

gnr2(lowa,1,ilowa_aut) = NaN;
gna2(lowa,1,1:nb) = gnaut2(lowa,1);
gna2(lowa,1,ilowa_rep) = NaN;

gnr2(hgha,1,ihgha_aut) = NaN;
gna2(hgha,1,1:nb) = gnaut2(hgha,1);
gna2(hgha,1,ihgha_rep) = NaN;

bpr2(lowa,1,ilowa_aut) = NaN;
bpr2(hgha,1,ihgha_aut) = NaN;

cnr2(lowa,1,ilowa_aut) = NaN;
cna2(lowa,1,1:nb) = cnaut2(lowa,1);
cna2(lowa,1,ilowa_rep) = NaN;

cnr2(hgha,1,ihgha_aut) = NaN;
cna2(hgha,1,1:nb) = cnaut2(hgha,1);
cna2(hgha,1,ihgha_rep) = NaN;

ctr2(lowa,1,ilowa_aut) = NaN;
cta2(lowa,1,1:nb) = ctaut2(lowa,1);
cta2(lowa,1,ilowa_rep) = NaN;

ctr2(hgha,1,ihgha_aut) = NaN;
cta2(hgha,1,1:nb) = ctaut2(hgha,1);
cta2(hgha,1,ihgha_rep) = NaN;

wwr2(lowa,1,ilowa_aut) = NaN;
wwa2(lowa,1,1:nb) = wwaut2(lowa,1);
wwa2(lowa,1,ilowa_rep) = NaN;

wwr2(hgha,1,ihgha_aut) = NaN;
wwa2(hgha,1,1:nb) = wwaut2(hgha,1);
wwa2(hgha,1,ihgha_rep) = NaN;

taur2(lowa,1,ilowa_aut) = NaN;
taua2(lowa,1,1:nb) = tauaut2(lowa,1);
taua2(lowa,1,ilowa_rep) = NaN;

taur2(hgha,1,ihgha_aut) = NaN;
taua2(hgha,1,1:nb) = tauaut2(hgha,1);
taua2(hgha,1,ihgha_rep) = NaN;

nettaxr2(lowa,1,ilowa_aut) = NaN;
nettaxa2(lowa,1,1:nb) = nettaxaut2(lowa,1);
nettaxa2(lowa,1,ilowa_rep) = NaN;

nettaxr2(hgha,1,ihgha_aut) = NaN;
nettaxa2(hgha,1,1:nb) = nettaxaut2(hgha,1);
nettaxa2(hgha,1,ihgha_rep) = NaN;

fiscalr2(lowa,1,ilowa_aut) = NaN;
fiscala2(lowa,1,1:nb) = fiscalaut2(lowa,1);
fiscala2(lowa,1,ilowa_rep) = NaN;

fiscalr2(hgha,1,ihgha_aut) = NaN;
fiscala2(hgha,1,1:nb) = fiscalaut2(hgha,1);
fiscala2(hgha,1,ihgha_rep) = NaN;

nfiscalr2(lowa,1,ilowa_aut) = NaN;
nfiscala2(lowa,1,1:nb) = nfiscalaut2(lowa,1);
nfiscala2(lowa,1,ilowa_rep) = NaN;

nfiscalr2(hgha,1,ihgha_aut) = NaN;
nfiscala2(hgha,1,1:nb) = nfiscalaut2(hgha,1);
nfiscala2(hgha,1,ihgha_rep) = NaN;

car2(lowa,1,ilowa_aut) = NaN;
car2(hgha,1,ihgha_aut) = NaN;

qr2(lowa,1,ilowa_aut) = NaN;
qr2(hgha,1,ihgha_aut) = NaN;

spr2(lowa,1,ilowa_aut) = NaN;
spr2(hgha,1,ihgha_aut) = NaN;

endb = 2;

nameFig = char([char(label_model) '_' char(label_graph) '_debt']);
figure('name',nameFig);
subplot(5,2,1);P = plot(b,100*squeeze(1-hr2(lowa,1,:)),b,100*squeeze(1-ha2(lowa,1,:)),b,100*squeeze(1-hr2(hgha,1,:)),b,100*squeeze(1-ha2(hgha,1,:)));
title('unemployment','Fontweight','normal');
AX1 = gca;box on;ylabel('%');
set(AX1,'YTick',0:10:30);
axis(AX1, [min(b) endb -1 30]);
set(P,NameArray,ValueArray);
subplot(5,2,2);P = plot(b,squeeze(pr2(lowa,1,:)),b,squeeze(pr2(hgha,1,:)),b,squeeze(pa2(lowa,1,:)),b,squeeze(pa2(hgha,1,:)));
title('p^{N}','Fontweight','normal'); 
legend('low y^T','high y^T','Location','NorthEast','Orientation','vertical');
AX1 = gca;box on;
set(AX1,'YTick',3.5:0.5:5.5);
axis(AX1, [min(b) endb 3.5 5.5]);
set(P,NameArray,ValueArray2);
subplot(5,2,3);P = plot(b,squeeze(gnr2(lowa,1,:)),b,squeeze(gna2(lowa,1,:)),b,squeeze(gnr2(hgha,1,:)),b,squeeze(gna2(hgha,1,:)));
title('g^{N}','Fontweight','normal'); 
AX1 = gca;box on;
set(AX1,'YTick',0.12:0.08:0.28);
axis(AX1, [min(b) endb 0.12 0.28]);
set(P,NameArray,ValueArray);
subplot(5,2,4);P = plot(b,-squeeze(bpr2(lowa,1,:))/(lambda+r),b,-squeeze(bpr2(hgha,1,:))/(lambda+r));
title('debt','Fontweight','normal');
AX1 = gca;box on;
set(AX1,'YTick',-1.6:0.8:1.6);
axis(AX1, [min(b) endb -1.6 1.6]);
set(P,NameArray0,ValueArray0);
subplot(5,2,5);P = plot(b,squeeze(cnr2(lowa,1,:)),b,squeeze(cna2(lowa,1,:)),b,squeeze(cnr2(hgha,1,:)),b,squeeze(cna2(hgha,1,:)));
title('c^{N}','Fontweight','normal');
AX1 = gca;box on;
set(AX1,'YTick',0.6:0.1:0.9);
axis(AX1, [min(b) endb 0.6 0.9]);
set(P,NameArray,ValueArray);
subplot(5,2,6);P = plot(b,squeeze(ctr2(lowa,1,:)),b,squeeze(cta2(lowa,1,:)),b,squeeze(ctr2(hgha,1,:)),b,squeeze(cta2(hgha,1,:)));
title('c^{T}','Fontweight','normal');
AX1 = gca;box on;
set(AX1,'YTick',0.8:0.2:1.4);
axis(AX1, [min(b) endb 0.8 1.4]);
set(P,NameArray,ValueArray);
subplot(5,2,7);P = plot(b,squeeze(wwr2(lowa,1,:)),b,squeeze(wwa2(lowa,1,:)),b,squeeze(wwr2(hgha,1,:)),b,squeeze(wwa2(hgha,1,:)));
title('wages','Fontweight','normal');xlim([min(b) max(b)]);
AX1 = gca;box on;
set(AX1,'YTick',3:0.5:4.5);
axis(AX1, [min(b) endb 3 4.5]);
set(P,NameArray,ValueArray);
%version 1: with bond spreads
subplot(5,2,8);P = plot(b,squeeze(spr2(lowa,1,:)),b,squeeze(spr2(hgha,1,:)));
title('spreads','Fontweight','normal');
set(P,NameArray0,ValueArray0);ylabel('%');
AX1 = gca;box on;
set(AX1,'YTick',0:20:40);
axis(AX1, [min(b) endb -2 40]);
%version 2: with current account
subplot(5,2,9);P = plot(b,-squeeze(car2(lowa,1,:)),b,-squeeze(car2(hgha,1,:)));
title('CA balance','Fontweight','normal');
AX1 = gca;box on;
set(AX1,'YTick',-0.3:0.3:0.3);
axis(AX1, [min(b) endb -0.3 0.3]);
set(P,NameArray0,ValueArray0);
%version 3: with (net) taxes
subplot(5,2,10);P = plot(b,squeeze(taur2(lowa,1,:)),b,squeeze(taua2(lowa,1,:)),b,squeeze(taur2(hgha,1,:)),b,squeeze(taua2(hgha,1,:)));
set(P,NameArray,ValueArray);
hold on;
PP = plot(b,squeeze(nettaxr2(lowa,1,:)),b,squeeze(nettaxa2(lowa,1,:)),b,squeeze(nettaxr2(hgha,1,:)),b,squeeze(nettaxa2(hgha,1,:)));
title('(net) taxes','Fontweight','normal');
AX1 = gca;box on;
set(AX1,'YTick',0.25:0.5:1.25);
axis(AX1, [min(b) endb 0.25 1.25]);
set(PP,NameArray,ValueArray5);
saveas(gcf, fullfile(plot_path, nameFig),format_chart);
end;
